Partial asynchrony of coniferous forest carbon sources and sinks at the intra-annual time scale

As major terrestrial carbon sinks, forests play an important role in mitigating climate change. The relationship between the seasonal uptake of carbon and its allocation to woody biomass remains poorly understood, leaving a significant gap in our capacity to predict carbon sequestration by forests. Here, we compare the intra-annual dynamics of carbon fluxes and wood formation across the Northern hemisphere, from carbon assimilation and the formation of non-structural carbon compounds to their incorporation in woody tissues. We show temporally coupled seasonal peaks of carbon assimilation (GPP) and wood cell differentiation, while the two processes are substantially decoupled during off-peak periods. Peaks of cambial activity occur substantially earlier compared to GPP, suggesting the buffer role of non-structural carbohydrates between the processes of carbon assimilation and allocation to wood. Our findings suggest that high-resolution seasonal data of ecosystem carbon fluxes, wood formation and the associated physiological processes may reduce uncertainties in carbon source-sink relationships at different spatial scales, from stand to ecosystem levels.

This study uses wood formation data collected in 81 sites belonging to boreal, temperate and Mediterranean biomes.All the data followed the criteria or procedures applied either in the field or laboratory, as described below and according to the methodology described by Rossi, Anfodillo, et al. (2006); Rossi, Deslauriers, et al. (2006).
The timings of wood formation were determined by monitoring healthy dominant trees.
The sample size ranged from 1 to 55 trees among all sites throughout the entire growing seasons of 1998 to 2018.Stem microcores were collected weekly, or occasionally biweekly, at breast height (i.e., 1.3 m) using surgical bone-sampling needles or a Trephor tool.The samples included mature and developing xylem of the current year, the cambial zone and adjacent phloem, and at least one previous complete tree ring.The microcores were fixed in propionic or acetic acid solutions mixed with formaldehyde and stored in ethanol-water at 5 °C.They were then dehydrated with successive immersions in ethanol and D-limonene and finally embedded in paraffin or glycol methacrylate (samples from Switzerland were not embedded).The microcores were cut with rotary or sledge microtomes to obtain 10 to 30 µm thickness transverse sections.
The sections obtained were stained with cresyl violet acetate or a mixture of safranin and astra/Alcian blue, then examined by light microscopy (bright-field and polarised light).
Cambial initials of the vascular cambium, a secondary meristem, divide outward and inward to produce phloem and xylem mother cells that, in turn, form new phloem and xylem tissues.The process of tracheid formation (i.e., cell differentiation) is divided into different phenological phases including cell enlargement, secondary cell wall thickening and lignification, and then programmed cell death, leading to the mature stage.Cambial cells are characterised by thin cell walls and small radial diameters.The enlargement zone is represented by the absence of glistening under polarised light, which indicates the presence of only primary cell walls.Cells undergoing secondary cell wall formation glistened under polarised light.In mature cells, Cresyl violet acetate reacts with lignin, turning from violet to blue.Maturation is reached when the cell walls are entirely blue.
To investigate the variations in climate among the locations where we monitored wood formation, we sourced bioclimatic data from the CHELSA bioclimatic database V2.1, offering a spatial resolution of 30 arcseconds.From the set of 19 bioclimatic parameters, we curated a selection of seven variables, as detailed in the Supplementary Table 1.We excluded parameters that delivered overly broad climate descriptions, like annual temperature and precipitation, and filtered out variables that exhibited strong autocorrelations (with a correlation coefficient exceeding |0.7|).To categorize our study sites based on their climate-related attributes, we applied the Partitioning Around Medoids clustering algorithm (PAM), an extension of k-means clustering approach.We employed the Within-Sum-of-Squares method (WSS) to identify the optimal number of clusters.This method minimizes the internal distances between data points within each cluster.We determined that four clusters best suited the task of grouping our wood formation study sites effectively (Supplementary Figure 1).We proceeded to perform a Principal Component Analysis (PCA) on the selected bioclimatic variables.This analytical technique allowed us to visualize and understand the climatic positioning of our 81 wood formation sites in relation to these variables.Pearson's correlation coefficient was employed to identify the climatic factors that influenced the ordering of wood formation study sites by principal components (Supplementary Table 1).

Bioclimatic variables
Abbreviation PC1 PC2 may not be directly comparable (Quentin et al., 2015).Ultimately, all data were normalized (min-max normalisation) on a scale from 0 to 1.

Flux data evergreen needleleaf forests (FLUXNET 2015)
We used tier-one level data from the FLUXNET2015 dataset (Pastorello et al., 2020) and extracted data at daily temporal aggregation from ENF (Evergreen Needleleaf Forests) sites.These sites consist of forest lands dominated by woody vegetation with a cover of >60% and height exceeding 2 meters.In addition, to reduce unwanted variability due to some specific characteristics of the site, we selected only data in which stands (1) belonged to boreal, temperate, or Mediterranean biomes, (2) were at least 15 years old, and (3) were not recently disturbed (e.g., burn sites).The dataset finally consisted of 39 sites.CO2 fluxes extracted for each site were Net Ecosystem Exchange (NEE), Ecosystem Respiration (RECO), and Gross Primary Production (GPP).Specifically, for the NEE we used the variable NEE_VUT_REF since it maintains the temporal variability and represents the ensemble (Pastorello et al., 2020).There are two main ways to estimate GPP and RECO (in units of g C m -2 d -1 ) by partitioning of the NEE computed with the variable Ustar (u * ) threshold (VUT): (1) night-time method (GPP_NT_VUT_REF and RECO_NT_VUT_REF); and (2) day-time method (GPP_DT_VUT_REF) (Pastorello et al., 2020).All results shown here use an average of the night-and day-time partitioning methods (RECO_NT_VUT_REF and RECO_DT_VUT_REF for RECO, GPP_NT_VUT_REF and GPP_DT_VUT_REF for GPP).Finally, C fluxes were grouped by biome, flux type (i.e., NEE, GPP, RECO), and site and normalized (min-max normalisation) on a scale from 0 to 1.

Flux data in wood formation sites (FLUXSAT V2.0)
Daily GPP data were extracted for each site and each year where the wood formation was monitored.Differently from the dataset obtained with the data from FLUXNET2015, in this case, while the samples were collected on coniferous species, the study area could consist of a mixed forest.These GPP products were extracted by FluxSat v2.0 (Joiner & Yoshida, 2021), where FluxSat refers to data derived using FLUXNET eddy covariance tower site data and the coincident satellite data.This dataset provides global gridded daily estimates of terrestrial GPP and uncertainties at 0.05° resolution for the period 01/03/2000 to the recent past (Joiner & Yoshida, 2021).The GPP was derived from the MODerate-resolution Imaging Spectroradiometer (MODIS) instruments on the NASA Terra and Aqua satellites using the Nadir Bidirectional Reflectance Distribution Function (BRDF)-Adjusted Reflectances (NBAR) product as input to neural networks that were used to globally upscale GPP estimated from selected FLUXNET 2015 eddy covariance tower sites.GPP value (in g C m -2 d -1 ) for each site resulted from the interpolation of GPP values of the four nearest pixels (Joiner & Yoshida, 2021).Supplementary Table 2 Alphabetical list of all the species monitored in the study.In total, 39 species of conifers belonging to 8 genera were the object of the study.

Supplementary Note 1: Bioclimatic Analysis
The PCA of the 81 sites in which wood formation assessments were carried out yielded eight components, of which two were significant, collectively accounting for 79.63% of the variance (Fig. S3 and Supplementary Table 1).1).
Partitioning Around Medoids clustering algorithm facilitated the distinction of four bioclimatic clusters within the Northern Hemisphere, with a clear border between climatic biomes (Fig. S3).The Mediterranean cluster showed direct correlations with the mean temperature of the driest, warmest, and coldest quarters (bio9, bio10, and bio11, respectively) while inversely correlating with precipitation seasonality (bio16) and precipitation of the driest and warmest quarters (bio17 and bio18, respectively) (Fig. S3).
The sites located in temperate biomes belonged to two main clusters, with the colder sites being directly correlated with precipitation seasonality (bio16) and precipitation of the driest and warmest quarters (bio17 and bio18, respectively), and inversely correlated with the mean temperature of the driest, warmest, and coldest quarters (bio9, bio10, bio11, respectively) (Fig. S3).The central temperate cluster was associated with most bioclimatic variables used in the analysis.Lastly, the boreal cluster was primarily determined by temperature seasonality (bio4) (Fig. S3).
biome, the peak is reached later than in other biomes, and the period of maximum activity (or minimum when considering soluble sugar concentrations) was shorter than in temperate and Mediterranean biomes that followed this specific order.Therefore, the amplitude of the oscillations among biomes for each fitting was larger in Mediterranean biomes than temperate and boreal biomes.
For each curve we estimated the timing of the maximum value (or minimum when considering soluble sugar concentrations) (Tables S4, S5, and S6).In addition to this, fitted curves pertaining of seasonal pattern of wood formation (i.e., cambial activity, cell enlargement, and cell wall thickening and lignification) and carbon fluxes derived from FluxNet data (i.e., NEE, GPP, and RECO) were utilized to estimate the timing of key percentiles, namely the 10 th , 25 th , 50 th , 75 th and 90 th percentiles, for both the ascending and descending portions of the curves.We opted to use the 75 th percentile as the threshold for defining the period of maximum activity (MA).However, when the curve exhibited a minimum in soluble sugar concentrations, the period of maximum activity was determined based on the 25 th percentile.
For each function, the area under the curve (AUC) was quantified, integrating the function itself (Tables S4, S5, and S6).The process of integration was also repeated by calculating the AUC of the period of MA for each function.This latter was obtained by using a definite integral with an interval consisting of a lower and upper limit represented by the points of intersection used to identify 75% of the maximum value of the curve and defining the period of MA (Tables S4, S5, and S6).Supplementary Figure 7 presents a comparison among the AUC of the period of MA for the processes of cell enlargement, GPP, and cell wall thickening and lignification.
Supplementary Table 3 Summary of all curve fitting.RSE indicates the residual standard error.All fittings were significant, at least with a p-value < 0.05.

Cambial activity and xylem formation
Cambial activity and wood formation stages were fitted with skewed normal distribution curves with a residual standard error (RSE) ranging between 0.14 and 0.26 (Supplementary Table 3).For cambial activity and cells differentiation stages, the peak was reached later in the boreal biome, and the period of maximum activity was shorter than in temperate and Mediterranean biomes that followed in this specific order (Supplementary Figure 4).Cambial activity in boreal and temperate ecosystems showed a peak at the beginning of summer and end of spring, respectively.In the Mediterranean biome, the peak was substantially earlier during spring.Cell differentiation showed a peak in the number of cells around the summer solstice (i.e., middle of June -start of July), with cell enlargement peaks localized in mid-June and cell wall thickening and lignification peaks localized at the beginning of July (Supplementary Figure 4).
The Mediterranean climate is characterised by dry and hot summers, and optimal growth conditions occur during the cooler and rainy springs and autumns.Being sensitive to drought, radial growth may cease in response to water shortage (Muller et al., 2011), leading to a bimodal growth pattern with a temporary cessation or reduction of radial growth in summer (Camarero et al., 2010).For this reason, we tried to fit a bimodal curve for the Mediterranean biome, but, finally, comparing the fittings, it was decided to use a skewed normal distribution that guaranteed a better fitting.Growth resumption in autumn after summer cease in response to drought is a plastic adaptation to large spatiotemporal variability in the climatic conditions in the Mediterranean region.However, even if the drought period can be temporally localized during a broad time window such as the summer season, the exact period when a moisture shortage comes into play can vary spatially and year by year in the Mediterranean region (Martins et al., 2012;Ruiz-Sinoga et al., 2012).Moreover, it is largely recognised that different species can show different responses to the dry period (Camarero et al., 2010).For all these reasons, in this context, by considering several species and different years, it was unlikely to highlight a bimodal growth pattern in the Mediterranean biome.
Supplementary Table 4 Maximum value (Max), onset, ending and duration of the period of maximum activity (MA), total and period of maximum activity area under the curve (AUC) and their ratio for all the curves describing the seasonal dynamics of wood formation (i.e., cambial activity, cell enlargement, and cell wall thickening and lignification).showed a residual standard error (RSE) ranging between 0.10 and 0.27 (Supplementary Table 3).The seasonal patterns of starch and soluble sugars were characterised by opposite temporal dynamics, particularly evident in needles with a very sharp starch peak in late spring-early summer strictly followed by the seasonal minimum for soluble sugars (Supplementary Figure 5).The amplitude of the oscillations among biomes was most prominent for starch and particularly noticeable in needles.Starch levels peaked belowground first (~early spring), then in stems (mid-spring) and finally in needles (late spring-early summer) but not regarding the boreal biome where the peak of starch in roots was later than the peak in needles and stem.Soluble sugars in roots and stems were less variable, with a hint of a seasonal minimum around late spring-early summer (Supplementary Figure 5).However, soluble sugar concentration in roots showed a peak in autumn in the Mediterranean biome.Boreal and temperate ecosystems showed contrasting temporal dynamics for starch and soluble sugar.Starch peaked in needles, stems and roots around late spring-early summer, mid-spring, and midsummer, respectively.In contrast, soluble sugars were lowest in all organs from late spring to midsummer.Temperate biomes were characterised by maximum starch concentrations towards late spring-early summer, particularly in needles, strictly followed by minimum levels of soluble sugars in all organs.In the Mediterranean, the magnitude of such oscillation was generally lower than those observed for boreal and temperate ecosystems (Supplementary Figure 5).

Supplementary Table 5
Maximum value (Max), onset, ending and duration of the period of maximum activity (MA), total and period of maximum activity area under the curve (AUC) and their ratio for all the curves describing the seasonal dynamics of NSC (i.e., starch and soluble sugars in needles, stem and roots).

Figure 2
Whittaker biome plot showing the comparison between the study site classification in the present study (i.e.boreal, temperate and Mediterranean biomes) and Whittaker biome classification according to mean annual temperature (°C) and annual precipitation (cm).